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ABSTRACT 

We investigate the evolution of globular clusters using A-body calculations and anisotropic 
Fokker-Planck calculations. The models include a mass spectrum, mass loss due to stellar 
evolution, and the tidal field of the parent galaxy. Recent A-body calculations have revealed 
a serious discrepancy between the results of A-body calculations and isotropic Fokker-Planck 
calculations. The main reason for the discrepancy is an oversimplified treatment of the tidal 
field employed in the isotropic Fokker-Planck models. In this paper we perform a series 
of calculations with anisotropic Fokker-Planck models with a better treatment of the tidal 
boundary and compare these with A-body calculations. The new tidal boundary condition 
in our Fokker-Planck model includes one free parameter. We find that a single value of this 
parameter gives satisfactory agreement between the A-body and Fokker-Planck models over a 
wide range of initial conditions. 

Using the improved Fokker-Planck model, we carry out an extensive survey of the evolution 
of globular clusters over a wide range of initial conditions varying the slope of the mass function, 
the central concentration, and the relaxation time. The evolution of clusters is followed up to 
the moment of core collapse or the disruption of the clusters in the tidal field of the parent 
galaxy. In general, our model clusters, calculated with the anisotropic Fokker-Planck model with 
the improved treatment for the tidal boundary, live longer than isotropic models. The difference 
in the lifetime between the isotropic and anisotropic models is particularly large when the effect 
of mass loss via stellar evolution is rather significant. On the other hand the difference is small 
for relaxation-dominated clusters which initially have steep mass functions and high central 
concentrations . 

Subject headings: Galaxy: kinematics and dynamics — galaxies: kinematics and dynamics — 
galaxies: star clusters — globular clusters: general — open clusters and associations: general — 
methods: numerical 
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1. Introduction 

The most reliable method to model the evolution of dense star clusters is the direct integration of the 
equations of motions of all stars via iV-body calculations. The high cost of such calculations, however, 
limits the choice of the number of stars to a few ten thousands, a small number compared to real globular 
clusters, and also limits the number of calculations which can be performed on a practical time scale. 

In the last few decades Fokkcr-Planck models have gained in popularity for studying the dynamics of 
globular clusters (see Meylan & Heggie 1997 for a recent review). The relatively low computational cost 
of Fokker-Planck models makes it possible to perform a large number of calculations over broad ranges of 
initial parameters. This open the possibility to study the initial parameter space which may have been 
available at the formation of the population of observed globular clusters. 

On the other hand, one has to pay for the higher speed of Fokker-Planck calculations with a less 
accuracy of these models. It is therefore important to examine the accuracy of Fokker-Planck models. One 
way of such examinations is to do a comparative study between Fokker-Planck and iV-body calculations. 

1.1. What happened before 

The development of approximate methods for modeling the dynamical evolution of star clusters started 
in the 1970's using Monte-Carlo methods (Spitzer 1987, and references therein). Several groups applied 
these methods to address various problems in stellar dynamics. Cohn (1979, 1980) developed a new 
method in which the orbit-averaged Fokker-Planck equation was solved directly using a finite-difference 
technique. Cohn (1980) demonstrated that his code could follow the evolution of the cluster into much 
more advanced stages of core collapse than Monte-Carlo codes (see, however, Giersz 1998, for a possible 
revival of Monte-Carlo methods). 

In Cohn's calculations globular clusters were modeled as fairly idealized systems, i.e.: isolated systems 
of identical point masses (Cohn 1979, 1980). This is quite a natural starting point for theoretical studies. 
After Cohn's seminal works, many contributions have been made to improve Cohn's Fokker-Planck scheme 
in terms of adding various astrophysical processes, such as a mass spectrum, primordial binaries, a tidal 
cut-off, gravitational shocks by the galactic disk and bulge, mass loss from stellar evolution, anisotropy of 
the velocity distribution, rotation of the clusters, etc. (see Meylan & Heggie 1997). 

An important landmark in modeling the evolution of globular clusters with Fokkcr-Planck calculations 
was set by Chernoff & Weinberg (1990, hereafter CW). Their models included a tidal cut-off by the parent 
galaxy, a stellar mass spectrum, and mass loss due to stellar evolution. 

CW boldly added simultaneously several new astrophysical phenomena to their models. However, 
generally speaking, each time a model is extended trying to make simulations more "realistic" , the number 
of assumptions usually increases; therefore thorough testing is required (not only for Fokker-Planck models) 
to eliminate the possible introduction of omissions and the violation of previously-made assumptions. It 
is not easy to assess the reliability of a model. The most reliable way to do this may well be detailed 
comparison with another model which is based on a completely different numerical method. 

Such comparisons have found good agreement between TV-body, Fokker-Planck, and gaseous models for 
isolated star clusters of point masses (Giersz and Heggie 1994a, 1994b; Giersz and Spurzem 1994; Spurzem 
and Takahashi 1995). These are important results because they demonstrate that two-body relaxation is a 
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main physical process which drives the dynamical evolution of star clusters, as was expected from theory 
(see, e.g., Spitzer 1987), and therefore that the time evolution of such relatively simple systems scales with 
the two-body relaxation time. 

Fukushige & Heggie (1995) and Portegies Zwart et al. (1998) demonstrated that iV-body calculations 
give significantly different results from the isotropic Fokker-Planck calculations of CW for some initial 
conditions. These iV-body and Fokker-Planck calculations included the steady tidal field of the parent 
galaxy and mass loss due to the evolution of single stars. The results of Portegies Zwart et al. (1998) 
can be summarized as follows: 1) The clusters simulated with the TV-body model live more than an order 
of magnitude longer than the comparable isotropic Fokker-Planck model; 2) The lifetime of the clusters 
depends on the number of stars N in a rather complex way. The clusters with a larger number of stars live 
shorter contrary to what was expected). A similar TV-dependence of iV-body models was also found in a 
collaborative experiment organized by Heggie (Heggie et al. 1999). 

Takahashi & Portegies Zwart (1998, TPZ hereafter) solved the discrepancies pointed out by Portegies 
Zwart et al. (1998) by using more general anisotropic Fokker-Planck models together with an improved 
treatment of the tidal boundary. The improved Fokker-Planck models are in in excellent agreement with 
iV-body models. 

The assumption of velocity isotropy and the oversimplified escape condition adopted by CW turned 
out to result in an enormous overestimate of the escape rate, i.e.: an underestimate of the life time of the 
cluster. TPZ took account of the orbital angular momentum as well as the orbital energy to determine 
escaping stars. TPZ also took account of that a star on an escape orbit requires some time to leave the 
cluster. The number of stars that escape per relaxation time depends therefore on the ratio of the relaxation 
time to the dynamical time, i.e.: on the number of stars. 

TPZ's implementation of the escape condition adds a free parameter v csc (denoted by a csc in TPZ) to 
Fokker-Planck models. This parameter determines the time scale on which escaping stars leave the cluster, 
and should be of the order of the crossing time at the tidal radius (Lee & Ostriker 1987; see also Ross 
et al. 1997). TPZ compared Fokker-Planck calculations with .V-body calculations to calibrated v csc for 
one specific set of initial conditions: density profile, mass function and relaxation time. However, it is not 
obvious that a single value of z/ osc gives good agreement between Fokker-Planck and .V-body models for a 
wide range of initial conditions. To test this hypothesis is the first main purpose of this paper. 

The second purpose of this paper is to carry out an extensive survey of the evolution of globular 
clusters in the Galaxy for various initial conditions. The same initial conditions will be used as those of 
the survey by CW. We perform the survey with the improved anisotropic Fokker-Planck models with the 
calibrated value for v csc . 

In the following section we review the anisotropic Fokker-Planck model and the .V-body model. In 
section || the stellar evolution models are described, and the initial conditions are reviewed in section 
|]. The results of a comparison between Fokker-Planck and TV-body models are presented in section 0. 
Fokker-Planck survey results are shown in section ||. Section ^ discusses our simulation results, including 
comparison with other iV-body simulations and with observations. We summarize our findings in section||. 



2. The Models 
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2.1. The Fokker-Planck model 

We use anisotropic orbit-averaged Fokker-Planck models (Cohn 1979; Takahashi 1995). The star 
cluster is assumed to be spherically symmetric and in dynamical equilibrium at any time. The distribution 
function / at time t is then a function of the energy of a star per unit mass, E, and the angular momentum 
per unit mass, J (see e.g., Spitzer 1987, section 1.2). On the other hand, isotropic Fokker-Planck models 
(Cohn 1980) assume that the distribution function does not depend on the angular momentum, i.e., that 
the velocity distribution is isotropic at every point in the cluster. 

Anisotropic Fokker-Planck models are more general than isotropic Fokker-Planck models and therefore 
preferable. The higher computational costs and more complex numerical implementation, however, caused 
anisotropic Fokker-Planck models to be less favored than their isotropic cousins (Cohn 1980). The numerical 
problem which was encountered by Cohn (1979) in his anisotropic Fokker-Planck calculations was solved 
by Takahashi (1995; 1996) mainly by implementing a new finite-difference scheme (see also Drukier et al. 
1999). 

All Fokker-Planck calculations in this paper are performed with the code which was developed by 
Takahashi (1995; 1996) and was later modified to include a mass function for single stars (Takahashi 1997), 
a tidal boundary (Takahashi et al. 1997; TPZ), and stellar evolution (TPZ). 

2.1.1. Determination of escaping stars 

TPZ demonstrated that the results of Fokker-Planck calculations are sensitive to the choice of the tidal 
boundary conditions. Two different criteria for the identification of escapers were tested by Takahashi et 
al. (1997) and by TPZ. These are: 1) the apocenter criterion, in which a star is removed if 

r a (E,J)>r t ; (1) 

and 2) the energy criterion, in which a star is removed if 

E> E t = -GM/r t . (2) 

Here r t is the tidal radius of the cluster, r a (E, J) is the apocenter distance of a star with energy E and 
angular momentum J, G is the gravitational constant, and M is the cluster mass. 

The apocenter criterion is considered more realistic than the energy criterion. The latter tends to 
remove too many stars too quickly, since it may even remove stars orbiting inside the tidal radius. For most 
calculations we adopt the apocenter criterion to determine escapers for Anisotropic Fokker-Planck models; 
we call these models Aa models. We also adopt the energy criterion for some calculations, these models are 
called Ae models. 



2.1.2. Introducing the crossing time in the Fokker-Planck model 

In the Fokker-Planck calculations performed by CW, escapers are removed instantaneously from the 
clusters. In the absence of stellar evolution this condition may be justified if the dynamical time scale of 
the cluster is negligible compared to the two-body relaxation time, i.e.; if the number of stars N — > oo. 
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TPZ demonstrated that the lifetime of stellar systems with small N is dramatically affected by the 
time scale on which escaping stars are removed from the cluster. Two treatments of escaping stars were 
tested by TPZ: 1) instantaneous escape: a star is removed from the cluster as soon as it satisfies an adopted 
escape criterion; 2) crossing-time escape: a star on an escape orbit hangs around in the cluster for a certain 
time before it really leaves the cluster. 

To implement the crossing-time escape condition in Fokker-Planck calculations we use the boundary 
condition of Lee & Ostriker (1987): 



dt ~ esc/ 



1 1/2 



(3) 



Here pt is the mean mass density within the tidal radius and v esc is a dimensionless constant which 
determines the time scale on which escapers leave the cluster. This equation is applied in a region of phase 
space where the adopted escape criterion is satisfied. Eq. (0) is derived assuming that escapers leave the 
cluster on a dynamical time scale. For TV — > oo with a fixed relaxation time, this condition becomes identical 
to removing escapers instantaneously, since the ratio of the crossing time to the relaxation time vanishes. 

When Eq. (|J) is used, the density at the tidal radius is generally still finite. The total mass of the 
cluster is then defined by the mass within r t . 



2.1.3. The mass function 

In our Fokker-Planck models a continuous mass function is represented by K discrete mass components. 
We use K = 20 components for all Fokker-Planck calculations. This number is considered to be sufficiently 
large for our calculations (see CW). At t — the mass of a star in mass-bin k is given by 

(m + \ (k -^ )/K 

m k = m_ — . (4) 

\m_ / 

Here m_ and m+ are the lower and upper limits of the initial mass function. The total mass in mass 
interval k is 

[ m K+i dN 

Mk — / — — mdm. (5) 
J dm 

Here dN / dm is the number of stars per unit mass interval. 



2.2. The TV-body model 

The iV-body portion of the simulations is carried out using the Kira integrator, operating within the 
Starlab (version 3.3) software environment (McMillan & Hut 1996; Portegies Zwart et al. 1998). Time 
integration of stellar orbits is accomplished using a fourth-order Hermite scheme (Makino & Aarseth 1992) 
Kira also incorporates block time steps (McMillan 1986a; 1986b; Makino 1991), special treatment of close 
two-body and multiple encounters of arbitrary complexity, and a robust treatment of stellar and binary 
evolution and stellar collisions (Portegies Zwart et al. 1999). The special-purpose GRAPE-4 (Makino 
ct al. 1997) system is used to accelerate the computation of gravitational forces between stars. The 
treatment of stellar mass loss is described in Portegies Zwart et al. (1998). A more complete description 
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of the Starlab environment is in preparation, but more information can be found at the following URL: 
tittp : //www. ias . edu/~ starlab. 



3. Stellar Evolution 
3.1. iV-body Models 

Stars in the TV-body calculations are evolved with the stellar evolution model SeBa^J (Portegies Zwart 
& Verbunt 1996, Section 2.1). In SeBa stars with mass m > 8M Q become neutron stars with mass 
m nn = 1.34 Mq, and lower mass stars become white dwarfs. The mass of these white dwarfs are given by 
the core mass of their progenitors. 

The evolution of the stars in the Fokker-Planck calculations which are compared with the TV-body runs 
is also computed using SeBa. In practice tabulated data as listed in Table |l| are used. 



3.2. Fokker-Planck models 

We use the same stellar evolution model as was adopted by CW when we compare our anisotropic 
Fokker-Planck models with CW's models. The stars in this stellar evolution model live somewhat shorter 
than those in SeBa (see Table [j]). 

Table | compares the stellar evolution model used by CW with SeBa. The main-sequence lifetimes 
in the CW model are obtained by fitting cubic splines to the data listed in Table 1 of CW. They do not 
mention the lifetimes of stars with mass m < 0.83Af Q . Some clusters become considerably older than 
15Gyr and it is necessary to use proper lifetimes for these stars as well. We assume that the lifetime of 
stars with m < 0.83M© is proportional to m -3 5 (Drukier 1995). This extrapolation gives also a good fit to 
the lifetime of the low-mass stars in SeBa (see Table [I]). 

At the end of the main-sequence life, a star with an initial mass m; n i forms a compact object with mass 
m nn . In the model adopted by CW the remnant mass is given by 

r 0.58 + 0.22(m ini - 1), for m ini < 4.7 [M©], 
m Rn =l 0, for 4.7 < m ini < 8.0 [M ], (6) 

[ 1.4, for 8.0 < m ini < 15.0 [M Q \. 

Stars with mass m; n i <J 0.46 Af@ do not lose any mass in a stellar wind as they turn into white dwarfs. 
The implementation in the Fokker-Planck model is as follows: The mass of the /c-th component (see 



2.1.3 ) is decreased from its initial mass m™ to the final (remnant) mass m^ n linearly with time t during 



the interval iMs(m^ 1 / 2 ) < t < iMs(m™ 1 ^ 2 ), where iMs(m) is the main-sequence lifetime of a star with 
mass m. 



3 The name SeBa is adopted from the Egyptian word for 'to teach', 'the door to knowledge' or '(multiple) star'. The exact 
meaning depends on the hieroglyphic spelling. 
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4. Initial Conditions 

The initial conditions for our Fokker- Planck and ./V-body calculations are identical to those of the 
survey performed by CW. It is assumed that the galaxy has a flat rotation curve and the star clusters are 
on circular orbits with v g — 220km/s around the galactic center. The tidal radius changes with time as 
rt(t) oc M 1 / 3 (i). The initial density profiles of the star clusters are set up from King models (King 1966). 
Calculations are performed with three different values of Wo = 1,3 and 7, where W is the dimensionless 
central potential of the King models. The initial mass function is given by 

dN oc m- a dm, (7) 

with a =1.5, 2.5 and 3.5 between to = 0.4 Mq and 15 Mq. 

In addition to the two dimensionless parameters W and a, there is one more free parameter: the 
initial relaxation time. Following CW, the mean relaxation time within the tidal radius is defined as 

(eq. [6] of CW with r = r t and c\ — 1), where to* is a representative stellar mass. This equation can be 
rewritten with setting to* = Mq as 

trix,cw = 2.57 x lOV^cw, (9) 

where 

_ M R s [220 kms- 1 ] 1 
FCW= [MoTIkFc] v e hJt (10) 

defines the cluster family (see CW). Here i? g is the distance to the Galactic center and v g is the circular 
speed of the cluster around the Galactic center (cf. eqs. [1] and [2] of CW). The crossing time of the cluster 
is then 

t aK ,(j±-\ cw £^®!l n jV. (11) 

\GM J llx ' ow M 1 ; 

A group of clusters of the same family (the same few) have the same initial relaxation time (^|). Our 
survey covers CW's families 1, 2, 3 and 4. Table | reviews the properties of the selected families. For 
convenience Table galso gives the more usual half-mass relaxation time (Spitzer 1987, p. 40), 

M 1 /2 r 3/2 

i rh = 0.138- h 



GV2( m ) In TV 



:i54xl0\vv7v,v(^) 3/2 ^, (12) 

7t ) (to) 



where is the half-mass radius and (to) is the mean stellar mass. 



5. Comparison between Fokker-Planck and TV-body Models 

We compare the results of anisotropic Fokker-Planck calculations (Aa models) and those of TV-body 
calculations. The escaping stars in the Fokker-Planck models are removed on a crossing time scale using 
Eq. (p|), but first we have to calibrate the value for u esc . 
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5.1. The calibration of v csc 

The rate of mass loss in the Fokker-Planck models changes monotonically with ^ csc ; larger v csc causes 
the models to lose mass on a shorter time scale. Increasing v osc has qualitatively the same effect on the 
mass loss speed as increasing N (i.e. increasing i r i x /t cr ) with is csc held fixed. 

The right choice of the initial conditions with which v csc should be calibrated is rather subtle (TPZ). 
We use a wide range of initial conditions to find the most appropriate value for v csc . No statistical analysis 
is performed for deciding which value of i> csc is the best, since the criterion with which a statistical test 
should comply is not clear. Rather we judge from global agreement in mass evolution by eye. 

Figure [j] shows the total mass evolution for Fokker-Planck calculations with v csc —2, 2.5, and 3 for 
selected initial models and compares them with iV-body calculations. These initial models were selected 
because stellar evolution as well as two-body relaxation play a significant role for these models. 

For the models shown in the top two panels, v csc — 3 gives good agreement between the Fokker-Planck 
and TV-body models at later epochs for N =1K; For N =16K and 32K, a value of z^ csc between 2 and 2.5 
is preferable, except for the final phase. For the models shown in the middle panels, v csc — 3 gives good 
agreement for the entire evolution. For the models shown in the bottom panels, the evolution is rather 
insensitive to v csc as well as N, and the agreement between the Fokker-Planck and iV-body models is good. 

From these comparisons we conclude that the best agreement between A^-body and Fokker-Planck 
models is achieved for 2 < z/ esc < 3, and we adopt the median value of u esc — 2.5 for the rest of our 
calculations. 



5.2. The time evolution of the total mass 

Here we present comparisons between A^-body and Fokker-Planck calculations for various initial 
conditions. A part of the results were already presented in the previous subsection, but are reproduced here 
with new information. 

TPZ performed Fokker-Planck calculations for the initial model of W — 3, a — 2.5 and family 1 with 
using i/ csc = 2. Fig. ||a shows repeated calculations with the new adopted value of u esc — 2.5. We also 
perform A^-body and Fokker-Planck calculations with initially 8K and 16K stars for other selected sets of 
initial conditions (for economical reasons one set of initial conditions is computed with only 4K and 8K 
stars); see Figs, ^b-f The Af-body calculations are continued up to 10 Gyr unless disruption occurs before 
this time. To see the convergence for N — > oo one extra Fokker-Planck calculation for each set of initial 
conditions is performed with the instantaneous escape condition. 

The dash-dotted lines in Fig. || represent the expected evolution of the total mass when no escapers are 
allowed, i.e., when mass loss occurs only through stellar evolution. They give upper limits to the cluster 
mass at each instance. (The A^-body model of Fig. ||a with IK stars seems to lose mass less quickly than the 
pure stellar-evolution mass loss case! But this is caused by random rendering in the small number of stars.) 

Fig. H shows the time evolution of the relative difference in the total mass between the Fokker-Planck 
and A^-body models, which we define as [Mp_p(t) — -^V-body WJ/'^A'-bodyW" The rapid increase of 
the relative difference seen in the upper two panels of Fig. || is the result of the disruption of the A^-body 
models. 
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Fig. [| plots the dimensionless mass loss rate £ as a function of time, where time is expressed in units of 
the initial half-mass relaxation time t^Q. The cluster mass loss rate dM/dt may be divided into two parts; 

dM _ /dM\ fdM" 
~dT ~ \ dTj esc + \dT 



(13) 



where the first term in the right-hand side represents the mass loss due to escapers from the tidal radius 
and the second term represents the mass loss due to stellar evolution. Since the stellar evolution is the same 
in the Fokker- Planck and TV-body models, it may be useful to separate the stellar evolution mass loss from 
the total mass loss in order to further clarify the difference between the two models. Thus in Fig. |] we plot 
the mass loss rate due to escapers £ osc as well as the total mass loss rate £ tot defined as follows: 



_^ho dM _trho f dM 

4tot " M dt ' tosc ~ M \ dt 



(14) 



When there is no stellar evolution, £ csc is usually much less than 1 (Lee & Goodman 1995). This is to be 
expected, since in that case mass loss proceeds only on two-body relaxation time scale. In the present 
models, since the stellar evolution mass loss causes the shrink of the tidal radius and thus produces more 
escapers, £ e sc exceeds 1 sometimes. 

Fig. H as well as Figs. || and |^ show that the Fokker-Planck and iV-body calculations agree fairly well 
over the entire range of initial conditions we investigated. Apparently the single value of v esc — 2.5 is 
applicable. The relative difference of the total mass at each time is generally less than 10%. However, for 
the cases of (Wo, a, family) = (3,2.5,4) and (7, 1.5,4), the difference is rather large. This large difference 
seems to well correlate with large £ osc (> 1), especially in the A-body models. Also when the cluster is 
about to dissolve (see panels [a] and [b]), £ eS c becomes larger than 1 and the Fokker-Planck and A-body 
models deviate. 

In the case of (Wo, a, family) = (3,2.5,4), the A-body models tend to lose mass more quickly than 
the Fokker-Planck models. The mass loss in the A-body models accelerates rapidly, compared with the 
Fokker-Planck models, from about 1 Gyr for the 16K model and from about 2 Gyr for the 8K model, and 
after some time the rate of mass loss decreases again. A similar behavior, with smaller amplitudes, is also 
observed for the other initial conditions. This phased mass loss is caused by a retardation effect in the 
A-body models. Each time the clusters lose mass by stellar evolution and by loosing stars from the tidal 
boundary it requires several crossing times to fill up the outer layers of the cluster and to recover dynamical 
equilibrium. In this period mass loss mainly comes from stellar evolution only, which results in a temporary 
decrease in the rate of mass loss. Since the crossing time is long (t CI ~ 1 Gyr) in the models shown in 
Fig. ^d, such a phased mass loss is clearly seen in this case. (In the Fokker-Planck models dynamical 
equilibrium is assumed.) The large values of £ csc result from the rapid stellar-evolution-driven mass loss 
and the long relaxation time of these models (£ esc measures the fraction of mass lost per relaxation time) . 

The models with (Wo, a, family) = (7, 1.5, 4) are characterized by a flat initial mass function and a high 
initial concentration. Fig. |^e shows that the cluster loses about 80% of its initial mass within the first billion 
years, and still remains bound. The mass loss is mainly driven by stellar evolution. The rather flat initial 
mass function causes about 33% (in number) of the stars to leave the main-sequence within 1 Gyr. The 
fraction of mass lost from the cluster due to stellar evolution alone is therefore about 62%. Nevertheless the 
cluster still survives beyond 10 Gyr. The cluster is saved by its initial high central concentration (W = 7). 
In the first 1 Gyr the halo of the cluster is stripped almost completely, but the compact core survives. 
Even in this extreme case the agreement between Fokker-Planck and A-body models is fair. The difference 
between the Fokker-Planck and A-body models generally builds up in the first billion years (when £ csc > 1), 
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and stays almost constant afterwards. The Fokker- Planck models assume that the potential change is 
adiabatic, but this assumption is clearly violated during the early rapid mass-loss period, since i cr ~ 1 Gyr. 

The main results of this subsection are summarized as follows: 1) Aa Fokker-Planck models with 
is csc — 2.5 agree fairly well with iV-body models in all the cases we investigated; 2) however, the agreement 
becomes worse when a large fraction of the total mass is lost on a crossing time scale. Such rapid mass loss 
conflicts with the assumption of the adiabatic potential change which the Fokker-Planck model is based on. 
In real globular clusters the number of stars is at least one order of magnitude larger than the number of 
particles used in the present A^-body simulations, and therefore we expect that the Fokker-Planck model 
should work well for them. 



5.3. The mass function 

At the start of each calculation the mass function is independent of the distance to the cluster center. 
Mass loss by stellar evolution and mass segregation due to the dynamical evolution of the star cluster cause 
the global mass function as well as the local mass function to change with time. After some time the mass 
function becomes noticeably a function of the distance to the cluster center. 

Figure || shows the mass function for the A-body and the Fokker-Planck calculations for the model 
with W = 7, a = 2.5 and family 1 at the age of 10 Gyr. At this moment the clusters has lost about 60% 
of its initial mass (see Fig. |2|c) . The differences between the global mass- functions of the A-body model 
and of the Fokker-Planck model are negligible. Even the partial mass functions for the main-sequence stars 
and for the compact objects (white dwarfs and neutron stars) agree excellently. At an age of 10 Gyr the 
mass function for the main-sequence stars is somewhat flatter than the initial mass function, but the shape 
of the mass function is still well represented by a single power-law. Also for some other initial models we 
compared the mass functions at t = 10 Gyr between Fokker-Planck and A-body models, and found good 
agreement in all cases. 

6. Fokker-Planck Survey Results 

In this section we present the results of our anisotropic Fokker-Planck calculations in which escaping 
stars are removed instantaneously from the cluster. These calculations give lower limits to the lifetimes of 
the clusters (see TPZ). We use the same stellar evolution model as CW in the survey simulations. 

Table | summarizes the survey results (Aa models) and compares them with the results obtained by 
CW (their Table 5). Following CW, the destiny of the cluster is classified as 'C for core collapse or 'D' 
for disruption. The mass of the cluster M cn d (in units of the initial mass) and the time t cn d (in Gyr) at 
the moment we stop the calculations (resulting in C or D) are listed in Table ||. The listed end time for 
collapsing clusters means the moment of core collapse. These clusters will survive beyond the core collapse, 
but their post-collapse evolution is not followed in our calculations. 

The disruption time as well as the mass at that instant are rather ill defined in Fokker-Planck 
calculations. By definition a cluster ceases to exist (disrupts) at the moment M = 0. However 
Fokker-Planck calculations usually break earlier. Mass loss often accelerates rapidly before the moment the 
cluster ceases to exist (see Fig. fy. In such a case the Fokker-Planck model fails to find the self-consistent 
solution of the Poisson equation, and breaks (see CW, § Vila). CW identified this moment as disruption 
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time. This may be associated with the moment when the cluster loses dynamical equilibrium (Fukushige & 
Heggie 1995). CW's definition, however, is not free from ambiguity particularly for the mass at disruption, 
since the above-mentioned mass-loss acceleration is related to the breaking of numerical calculations. 

For clusters that dissolve we continued the calculations until the mass evolution curves start to bend 
downwards very rapidly (see Figure |^) . We define i V ert as the moment when the mass loss rate becomes 
virtually infinite. Then, however, the mass at i V ert is not well defined. Therefore, for convenience, we 
decided to take i cn d = 0.99i V ort as the disruption time, and then M cn( \ = M(t e nd)- We confirmed by eye 
that the point (t e nd, -^end) determined by the above definition corresponds to the "turn-off point" of the 
mass evolution curve. 



6.1. The total mass and the concentration as a function of time 

Figure || shows the evolution of the total mass for all the initial models. The circles plotted at the ends 
of some of the lines indicate the moments when those clusters experience core collapse. The destiny of the 
clusters is determined by the initial conditions (W , a, and Few)- The early evolution is strongly affected 
by the choice of W and a, and in less extent by the initial relaxation time, simply because the latter 
requires more time to affect the dynamics of the stellar system. The lifetimes of the clusters which are 
disrupted within ~ 1 Gyr, for example, are almost independent of the initial relaxation time. The evolution 
of these clusters is mainly driven by stellar mass loss. This causes the clusters to expand after which the 
Galaxy gobbles up the clusters. Core collapse occurs only in some long-lived (> 3 Gyr) clusters. For these 
clusters the mass loss rate gradually decreases with time (see also Portegies Zwart et al. 1998). 

Figure |?] shows the evolution of the concentration parameter c, which we define by (cf. King 1966) 

c = logio— • (15) 



Here the core radius r c is 



( 3<p \ 1/2 
WGpoJ 



(16) 



where po is the central density and v m o is the density-weighted velocity dispersion in the cluster center. 
Note that this core radius r c and the King core radius (King 1966) are not identical even for single-mass 
systems, since the central velocity dispersion is not used in the definition of the King radius. The two 
quantities get closer for more concentrated (higher W Q ) clusters (cf. Binney & Tremaine 1987, p. 235). 

In some cases (for W — 7, a = 1.5 of all the families, and for W = 1, a = 3.5 of family 1) the cluster 
experiences core collapse after almost all (~ 99%) of the initial mass is lost. We identify these clusters as 
collapsing clusters, because the concentration c rapidly increases at the final epochs (see Fig. [?]). If the 
initial mass of a cluster is relatively small (e.g. Mq <^ 10 4 Mq), the definite evidence of core collapse will be 
very hard to observe, since the number of stars in the core will become too small (e.g. N COIC <J 10). Deep 
core collapse in such systems can easily be prevented by the formation of binaries by three-body processes 
or by burning primordial binaries. Such processes are more effective in small N clusters. 

For W = 7, a = 1.5, irrespective of the families, the clusters have very small masses (compared to 
their initial masses) around t =10 Gyr. This implies that they may be observed at present as tiny globular 



clusters or as old open clusters (see § |7.4|) 



For all the models with a = 1.5, mass loss due to stellar evolution causes the concentration to decrease 
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with time before core collapse occurs. For a —2.5 and 3.5, the concentration does not decrease so much 
but stays almost constant. If, at any time, the concentration becomes smaller than c ~ 0.4 (compare with 
c = 0.58 for the Wo = 1 King model), the cluster quickly dissolves. As CW pointed out (see their section 
VIIc), the (scaled) density profiles just before disruption are always similar regardless of initial models. 

6.2. Comparison with Chernoff &: Weinberg (1990) 

6.2.1. Clusters that disrupt 

All our clusters of model Aa live longer than those of CW for the same initial conditions. The largest 
(relative) discrepancy is about a factor of ten which is found for (W a ,a) = (1,2.5) and (3,2.5). For the 
models with (W ,a) = (1, 1.5) and (3, 1.5), the discrepancy is about a factor of two. In the latter models 
the mass functions are very flat and the potentials are shallow, and thus stellar mass loss alone is enough to 
disrupt these clusters. It is therefore not surprising that even the Aa models cannot live so long compared 
with the CW models. The disruption time hardly depends on the initial relaxation time for those clusters 
which disrupt within ~ 1% of the half-mass relaxation time. 

Table |] compares our Ae and Aa models with CW's models for family 1. In addition to the data 
already given in Table |^, i en d and M cna \, two new quantities, t(M cn ^ cw) and M (i C nd,cw)> are presented. 
Here M en( j cw an( i £end,cw denote the final mass and time, respectively, of each of CW's models which are 
given in Table |^. These two quantities are particularly useful for comparing different models for disrupted 
clusters. 

The energy criterion for the tidal boundary is used for the CW models as well as for the Ae models. 
The Ae models are in this sense in between the CW and Aa models. Comparison with the Ae models helps 
us to understand how the apocenter criterion affects the results. The Ae models lose mass faster than the 
Aa models in all the cases shown in Table |], simply because the energy criterion introduces a larger escape 
region in phase space than the apocenter criterion (Takahashi et al. 1997; TPZ). The difference between 
the CW and Ae models is rather small in all cases. For each set of initial conditions the two models have 
the same destiny, disruption or collapse. The difference in the disruption time is within a factor of two, and 
that in f (-M en d,cw) is even smaller. We therefore conclude that the introduction of the apocenter criterion 
plays the most important role in extending the cluster lifetime. 

For (W , a) = (1, 2.5) and (3, 2.5) we find large differences between the CW and Aa models not only in 
t cn d and -M on( j but also in t(M e nd,Cw) and £(t e nd,cw)- This assures that the discrepancy between the two 
models are not merely due to differences in the definitions of the disruption time and mass. The clusters 
survive beyond the disruption times given by CW, which is also supported by A-body simulations (TPZ; § 
1- 

6.2.2. Clusters that collapse 

About half (seventeen) of the clusters calculated with Aa models experience core collapse, while in the 
calculations of CW less than 30% (ten) of the clusters experience core collapse. This means that seven of 
the models of CW are disrupted while the corresponding Aa models experience core collapse. In these Aa 
models core collapse occurs much later than the disruption time of the CW models; a too high rate of mass 
loss in the CW models disrupts the clusters too quickly. 
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All the clusters which experience core collapse according to CW do experience core collapse according 
to our calculations also. Core collapse occurs after two-body relaxation starts to dominate the evolution of 
the clusters, i.e.; when stellar-evolution mass loss has sufficiently decreased. Since CW's models generally 
overestimate the rate of mass loss, all models that finally develop a collapsed core according to their 
calculations experience core collapse in our calculations too. 

For (W , a) = (7,2.5) and (7,3.5) of all the families, both the models of CW and our Aa models 
experience core collapse. For these models, the moment of core collapse approximately scales with the 
initial relaxation time. This indicates that their evolution is dominated by two-body relaxation, not by 
stellar evolution. In those cases, the differences in the moment of core collapse between CW and Aa models 
are less than 10 %. 

The introduction of anisotropy into the models itself does not have very significant effects on the 
moment of core collapse (see Takahashi 1995; Takahashi et al. 1997). Core collapse proceeds almost 
independently of the outer parts of the cluster and the anisotropy in the core is always small. The difference 
between the energy and apocenter criteria affects mass loss from the tidal boundary, and the mass loss 
affects core collapse time by reducing the relaxation time. But the effect of the mass loss on core collapse 
time is only significant for low-concentration clusters which lose a substantial fraction of their mass before 
collapse (cf. Quinlan 1996). 

In addition, when mass loss proceeds slowly on the relaxation time scale, the difference between the 
apocenter criterion and the energy criterion is less important than when mass loss proceeds more rapidly 
driven by stellar evolution. The stars which exit in a region lying between the boundary of the energy 
criterion and that of the apocenter criterion in phase space (see, Takahashi et al. 1997, Fig. 1; or TPZ, 
Fig. 1) are responsible for slower mass loss in apocenter-criterion models. These stars diffuse over the 
tidal boundary on a time scale of the order of the relaxation time. If mass loss by stellar evolution is 
ineffective and the cluster evolves slowly on the relaxation time scale the difference between apocenter- 
and energy-criterion models will be small (cf. Takahashi et al. 1997). On the other hand, if the cluster 
loses mass rapidly due to stellar evolution compared with the relaxation time scale, the above-mentioned 
phase-space region will play an important role as a mass reservoir. 

We therefore conclude that the moment of core collapse is correctly calculated with isotropic 
Fokker-Planck models as well as with anisotropic Fokker-Planck models as long as the cluster is 
relaxation-dominated. 

Figure ^ gives a graphical representation of Table ^ for family 1. The largest discrepancies occur around 
a boundary which separates disrupted clusters and collapsing clusters; the disrupted clusters are on the 
lower left region and the collapsing clusters are on the upper right region in this figure. 

For Wq = 3 and a = 3.5 of family 1 and 2, the differences in the collapse time between isotropic and 
anisotropic models are rather large. These models are near the separating boundary mentioned above; the 
CW models of family 1 and 2 collapse but those of family 3 and 4 dissolve. In such cases, even a small 
difference in the mass loss rate may cause a large difference in the fate of the cluster. 

6.3. Fokker-Planck Survey Results for N = 3 x 10 5 

Real star clusters contain a finite number of stars, typically N ~ 10 5 - 10 6 . It is not clear whether this 
number is large enough to be treated as infinite, though it is usually assumed in Fokker-Planck calculations. 
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For example, it is not clear whether the instantaneous escape condition, which is valid for N — > oo and is 
adopted in our survey simulations, is a good approximation for typical globular clusters. To test the effect 
of a finite crossing-time, we performed a series of calculations by using the crossing-time removal condition, 
Eq. (|), with N = 3 x 10 5 . 

Figure [)] presents the mass evolution for models of family 4 with N = 3 x 10 5 initially. Family 4 has 
the largest relaxation time and thus has the largest crossing time for fixed N. The initial total masses of 
these models are 7.3 x 10 5 M Q , 3.0 x 10 5 M©, and 2.0 x 10 5 M Q , and the crossing times t CI are 0.03 Gyr, 
0.06 Gyr, and 0.10 Gyr, for a — 1.5, 2.5, and 3.5, respectively. These models are compared with the models 
with the instantaneous escape condition (see Fig. ^). 

The relative difference between the models with JV = 3x 10 5 and those with N — > co for the clusters 
which dissolve within a few Gyr is rather large. In these cases the models with N = 3 x 10 5 live a factor of 
two or three longer than the N — > oo models. The lifetime of these clusters is comparable to the crossing 
time. The evolution of these models is therefore very sensitive to the way escapers are removed from the 
clusters; crossing-time removal or instantaneous removal. For long-lived clusters that survive for a Hubble 
time the details of how escaping stars are removed become less important; the instantaneous removal 
condition gives a good approximation for these cases. 



7. Discussion 

7.1. The anisotropic Fokker- Planck model 

We have shown that our anisotropic Fokker-Planck (Aa) models are in excellent agreement with 
A-body models for a wide range of initial conditions. Although expected from theoretical considerations, 
it is still quite striking that a single value of v esc suffices to achieve agreement between the very different 
types of numerical models and over such a wide range of initial conditions. Subtle differences between the 
Fokker-Planck and A-body results, however, remain. These are probably due to the greater detail and 
higher accuracy of the A-body models compared to the Fokker-Planck models. (The absence of binary 
formation via three-body processes in the Fokker-Planck models may also play a role here.) In some cases, 
the relatively simplistic implementation of the tidal field via a cut-off radius causes a somewhat peculiar 



behavior in the iV-body models (see § 7.2). Fokker-Planck models lose their credibility when the time scale 
for mass loss becomes comparable to the dynamical time scale, since in that case the assumption of an 
adiabatic change in the potential is violated. 

Concerning the models shown in the present paper, all the Aa models live longer than the isotropic 
models. However, there are cases where anisotropic models live longer than isotropic models. For example, 
Takahashi et al. (1997) show the results that, in the absence of stellar evolution, Aa and Ae models 
evaporate faster than an isotropic model in the tidal field. The shorter lifetime of these anisotropic models 
is induced by two-body relaxation at the inner regions which causes a quick emergence of elongated-orbit 
stars in the halo. The apocenters of these stars are frequently outside the tidal radius. Especially when 
the tidal field is relatively weak mass loss in anisotropic models is quicker than in isotropic models (see 
Takahashi & Lee 1999 for details). 

The apocenter criterion is considered more realistic than the energy criterion, but also this approach 
has its limitation. The implementation of a spherical cutoff radius instead of a self-consistent tidal field, 
for example, limits the conditions on which escapers are identified. However, also in a real cluster a 
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star must pass the tidal radius in order to leave the cluster. The apocenter criterion is therefore an 
important improvement over the simple energy criterion used in previous Fokker-Planck calculations. The 
implementation of the apocenter criterion, however, requires anisotropic Fokker-Planck models; isotropic 
Fokker-Planck models do not contain sufficient information about the stars' orbits to apply the apocenter 
criterion. Unfortunately it seems hard to adjust the energy criterion used in isotropic models such that the 
mass loss rate becomes similar to that for the apocenter criterion. This is because the apocenter criterion 
causes the cluster's mass-loss rate to depend on the angular-momentum distribution (anisotropy) of the 
stars near the tidal boundary, and the anisotropy behaves in a rather complex way: it changes with time, 
and depends on the mass of the stars and the distance to the cluster center (Takahashi 1997; Takahashi et 
al. 1997). 



7.2. Effect of the tidal field 

To investigate the effect of the self-consistent tidal field compared with the adopted tidal cut-off, we 
perform several A-body calculations with different implementations of the tidal field. The same initial 
conditions as for the models in Figure ||a are adopted. For each model, calculations are performed ten times 
using different sets of random numbers with 1024 stars. The results are presented in Figure [l(]. 

The solid line in Figure [w] represents the mean of ten calculations with a simple tidal cut-off (model 
TCR1, the c/2 deviations from the mean are shown in Figure ||a). The other lines represent the results 
for the calculations with a self-consistent tidal field. For all the tidal-field models the same tidal force is 
used but different cut-off (maximum) radii are selected. Stars in the tidal-field models are removed as soon 
as they go beyond r t (model TFR1; the T stands for tidal, F for field and the R stands for the radius, in 
units of the tidal radius, at which stars are removed from the simulation), 2r t (model TFR2), 10 r t (model 
TFR10), or 100 rt (model TFR100) from the cluster center. For computational reasons it is convenient to 
set a cut-off radius even for the calculations where a tidal field is used. The mass given in Figure is the 
mass contained in stars which are bound to the cluster. 

All the models presented in Figure |l0| are indistinguishable for the first one billion years, where mass 
loss is dominated by stellar evolution. After that, model TCR1 and model TFR1 start to lose mass more 
rapidly than the other models. This acceleration in mass loss is related to a too sudden removal of escaping 
stars which is caused by the small cut-off radius. Model TFR1 keeps losing mass at a rate slightly higher 
than the other tidal-field models. After about 2 Gyr the mass loss rate of model TCR1 becomes smaller 
than for the models with a tidal field. Not surprisingly, stars are removed more quickly if a tidal field is 
incorporated and stars are removed at r t . Figure [w] shows that an artificial cut-off radius in excess of 2r t 
hardly affects the results in the models where a self-consistent tidal field is incorporated. 

The Fokker-Planck model with IK stars (see Figure ^a) loses more mass in the first billion years than 
the A-body models (TCR1, TFR1 to TFR100). This difference may be caused by the assumed adiabatic 
response of the Fokker-Planck models to changes in the potential (Cohn 1979). The cluster loses about 
30% of its initial mass within 1 Gyr. Since the initial half-mass crossing time is about 0.5 Gyr, this change 
in the potential is quite impulsive. In spite of this breakdown of one of the fundamental assumptions in 
Fokker-Planck models, the agreement with the A-body model is, except model TFR1, fairly good. 

We have demonstrated that A-body and Fokker-Planck models with the tidal cut-off behave in a 
similar fashion as ./V-body models with the self-consistent tidal field, for one selected set of initial conditions. 
However, it is not obvious whether tidal cut-off models are always similar to tidal-field models. Giersz & 
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Heggie (1997) also investigated the difference between these two kinds of TV-body models. They concluded 
that the mass evolution is not much affected by the difference between the two implementations of the tidal 
boundary conditions. We like to encourage further investigations on this issue. 



7.3. Comparison with similar TV-body surveys by other authors 

Recently Aarseth & Heggie (1998) investigated the time scaling of TV-body models with respect to the 
number of stars. They incorporated the Galactic tidal field in the same way as Fukushige & Heggie (1995) 
did. The same stellar evolution tracks as used by CW were adopted. 

Aarseth & Heggie (1998) introduced "variable time scaling" which couples the stellar evolution time 
scale to the time scale on which the star cluster evolves dynamically. During the early stages of the 
evolution of the cluster, when the effect of stellar evolution is dominant, time scaling based on the crossing 
time is applied. At later stages, when two-body relaxation becomes dominant, the relaxation time is used 
for time scaling. This variable time scaling has a qualitative physical basis, but there remains uncertainty 
in when and how to change from one scaling to another. 

Using the variable time scaling Aarseth & Heggie (1998) find good agreement between their TV-body 
calculations and the Fokker-Planck calculations of CW for those clusters that survive long enough for their 
evolution to become dominated by two-body relaxation. For quickly disrupted clusters, the models of CW 
live too short compared to the models of Aarseth & Heggie (1998). These conclusions are qualitatively 
consistent with our conclusions on the comparison between our Fokker-Planck models and the CW models. 
A detailed comparison between Aarseth & Heggie's (1998) TV-body models and our Fokker-Planck models 
is complicated by the difference in the implementation of the tidal field. 

In Table | our Aa models are compared with the TV-body results of Aarseth & Heggie (1998, Table 
3) as well as those of Fukushige & Heggie (1995, Table 3) for family 1, in terms of core collapse time and 
mass, and disruption time and mass. The models of the former authors are denoted by "AH" , and those 
of the latter by "FH" . A number associated with each of these model names indicates the galactocentric 
distance R g (in kpc) at which that model is placed. The initial cluster mass is 5.45 x 10 4 M© at ~ lOkpc, 
and 1.49 x 10 5 M Q at ~ 4 kpc. The Aa models use the instantaneous escape condition and their evolution 
is, within one family, independent of R g Fukushige & Heggie (1995) use fixed time scaling based on the 
crossing time, which is considered to be appropriate for the early evolution of clusters. Table shows only 
those models for which they obtain definite lifetimes. 

The disruption time and mass, idis and Mdi s , are listed in Table |^. The TV-body models disrupt at the 
moment the mass has dropped to zero, but, as we discussed above, the Fokker-Planck calculations usually 
stop well before the mass goes to zero. Therefore we should be aware that idis for the iV-body models and 
that for Fokker-Planck models do not have exactly the same meaning. 

For (Wo, a) — (1, 2.5) and (3, 2.5) the results of CW are very different from those of FH and AH, but 
our Aa models are consistent with the results of FH and AH. For models with the flat initial mass function 
of (Ho, a) — (1, 1.5) and (3, 1.5), the results of CW deviate largely from the calculations of FH and AH. For 
these cases, the lifetimes of our Aa models (^0.02 Gyr) are about twice as long as those of CW's models, 
but they are still shorter than the lifetimes of the models of FH and AH by a factor of a few. This difference 
will be reduced when the crossing-time removal condition is applied for Fokker-Planck models. As discussed 



in § 3.3, for such short-lived clusters with small TV, the instantaneous removal condition gives somewhat 
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too short lifetimes (see Fig. 

The Aa models (as well as CW's models) with (Wo, a) = (7,2.5) and (7,3.5) experience core collapse 
at about the same time and at about the same cluster mass as the models of AH. For (Wo, a) = (1,3.5) 
and (7, 1.5) the Aa models collapse while the TV-body models do not. The masses of these Aa models at the 
moment of core collapse, however, are only ~ O.OlMo, and thus it is not surprising that the core collapse 
does not occur in the ./V-body systems comprising initially of only a few thousand particles. 



7.4. Comparison with the observations 

For a better comparison, all our calculations are performed with the same initial conditions as adopted 
by CW. These initial conditions, however, are probably not the best choice for real globular clusters. 
Recent observations have revealed that mass functions are ill represented by a single power-law (e.g., 
Scalo 1986; Kroupa et al. 1990) and that the lower stellar mass limit should be <J 0.1M® (e.g, De Marchi 
& Paresce 1995a; 1995b; Marconi et al. 1998; see also Chernoff 1993). Another uncertainty, which may 
seriously affect the simulation results, is the initial compactness of star clusters. The King models from 
which wc selected the initial density profiles may not be the best choice for the initial density profiles 
for tidally-truncated star clusters (Heggie & Ramamani 1995). It is not even clear whether real globular 
clusters were born filling up their tidal lobes. Initial under-filling of the tidal lobe may extend the cluster 
lifetime considerably Regardless of these uncertainties in the initial conditions, it will be still worthwhile 
to discuss some implications for the observed population of globular clusters. 

Our calculations indicate that clusters which are observed in a state of core collapse at present must 
have been born fairly concentrated, otherwise core collapse would not be occurring within the age of the 
Universe. At core collapse most of the models have lost £ 50% of their initial mass. This indicates that 
observed globular clusters in a state of core collapse were born with a high concentration and considerably 
more massive than today's mass. 

All the clusters with a concentration parameter c « 1 at t « 10 Gyr appear to be in a local minimum 
in concentration. These clusters recover their concentrations later as core collapse follows. A cluster which 
reaches a concentration of c w 0.4 at any time evaporates quickly thereafter. The least concentrated 
globular clusters in our Galaxy have a King concentration parameter of CKmg ~ 0.5 (Tragcr et al. 1993). 
This is consistent with our findings 

According to our survey most of globular clusters with a rather flat initial mass function (a = 1.5) 
disrupt within a Hubble time. This indicates that the majority of the present-day clusters probably had 
steeper initial mass functions. However, we cannot exclude the possibility that some of the globular 
clusters were born with an initial mass function of a ~ 1.5 but have still survived owing to the initial high 
concentration. For example, the model with W = 7, a = 1.5 and family 3 loses about 99% of its initial 
mass and has still a high concentration at 12 Gyr (see Figs. ^ and 0). The old open cluster Berkeley 17 has 
similar characteristics as this model, as we show below. 

Phelps (1997) derives an age for Berkeley 17 of 10 Gyr to 13 Gyr (but Carraro et al., 1999, estimates 
an age of 9 ± 1 Gyr), a distance to the Galactic center of i? g = 11 kpc, a size of ~ 5 pc, and a total mass in 



4 Note, however, that cxing is n °t identical to the c (see § and that the one-tr error in the observed cxing i s about 0.2 
(Trager et al. 1993). 
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visible stars of about 400 M fl 



The core of the model with W = 7, a = 1.5 and family 3 collapses at an age of 12 Gyr. At this moment 
the total mass is only 1.3% of its initial mass (see Table ||). Substitution of i? g = 11 kpc and v g — 220 km/s 



into Eq. (10) results in an initial mass of Mo ~ 2.3 x 10 5 M©. At the moment of core collapse the cluster 
has then a mass of ~ 3000 M Q and the half-mass radius of ~ 4 pc. The majority of this mass is hidden in 
stellar remnants, i.e.: ~ 9OOAf in neutron stars and ~ 1700M Q in white dwarfs, and only ~ 400M Q is in 
main-sequence stars (see also Heggie & Hut 1996; Vesperini & Heggie 1997 for dark remnants in globular 
clusters). This model cluster looks similar to Berkeley 17 in the size and luminous mass. 

We note that the fraction of neutron stars in our models is overestimated, since many are likely to 
escape from the cluster by a "kick" at the time of a supernova explosion. An intrinsic asymmetry in the 
supernova is believed to give neutron stars a kick velocity which is generally much higher than the escape 
velocity from the cluster (Drukier 1996; Portegies Zwart & van den Heuvel 1999). This effect is not taken 
into account in our calculations. 

Figure |ll| shows the global mass function for the model with W = 7, a = 1.5, and family 3, at 
t = 12 Gyr. The filled circles and the open circles represent the mass functions for the white dwarfs and for 
the main-sequence stars, respectively. The mass function for the main-sequence stars peaks at m~ 0.8Mq. 
The stars with m £ 0.8M@ have evolved off the main-sequence. The peculiar shape of the main-sequence 
mass function for m £ 0.8M Q , which decreases with decreasing m, is the result of the preferential loss of 
low mass stars through the tidal boundary; high mass stars are more centrally concentrated than low mass 
stars due to mass segregation and therefore have less chance to escape. In general the global mass function 
becomes flatter as the cluster loses mass (e.g. Vesperini & Heggie 1997). However, it is rare to find such 
an inverted global mass function as is shown in Fig. |ll| in our surveyed models. Such a mass function 
only appears in clusters which have lost ^ 99% of their initial mass on the relaxation time scale. These 
"remnant" clusters were formerly the cores of more massive clusters. High mass stars are over-abundant in 
such cores as a result of mass segregation. 

The mass function of the globular cluster NGC 6712 around the half-light radius also peaks near 
0.75 Mq (De Marchi et al. 1999), and looks similar to the mass function shown in Fig. [H]. (Although Fig. |ll| 
shows the global mass function, the local mass function at the half-mass radius is similar to that.) 

De Marchi et al. (1999) interpreted the observed inverted mass function as evidence for severe tidal 
disruption, which is consistent with the above model. Thus we suspect that NGC 6712 was born far more 
massive than what it is today. With the current mass of ~ 10 5 M Q (Pryor & Meylan 1993), its initial mass 
would be ~ 10 M@, if that model is applied. However, this prediction of the initial mass is rather uncertain 
quantitatively, because our survey is limited in many respects. The range of initial conditions covered by 
our survey is still limited, as mentioned above. Further we assume a steady tidal field. The orbit of NGC 
6712 is, in fact, likely to be very eccentric; the cluster is presently at i? g = 3.5 kpc, but the minimum and 
maximum galactocentric distances are predicted to be about 0.3 kpc and 7 kpc, respectively (Daupholc et 
al. 1996). Therefore, bulge and disk shocks (see Gnedin & Ostriker 1997; Murali & Weinberg 1997; and 
references therein) have probably affected the mass evolution of this cluster considerably, as pointed out by 
De Marchi et al. (1999). 

Since we argued above that Berkeley 17 may be represented by the model shown in Fig. [ll], it will be 
interesting to observe the mass function of Berkeley 17 down to very low mass stars. 
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8. Conclusions 

We have investigated the evolution of star clusters including the effects of the tidal field of the parent 
galaxy and mass loss from stellar evolution by using anisotropic Fokkcr-Planck models as well as TV-body 
models. Our new Fokker-Planck models, which include a new implementation of the tidal field, agree much 
better with iV-body models than isotropic Fokker-Planck models. 

The new implementation of the tidal boundary consists of two parts: the apocenter criterion assures 
that only those stars are removed from the stellar system for which apocenter radius exceeds the tidal 
radius of the cluster; the crossing-time removal condition prevents escaping stars from being removed in an 
infinitesimal time step, but keeps them in the stellar system for some time of the order of the crossing-time. 
The crossing-time removal introduces an additional free parameter, v csci to the Fokkcr-Planck model. 
We find that a single value of i^ csc = 2.5 provides good agreement between iV-body and Fokker-Planck 
calculations over a wide range of initial conditions. These improvements to the Fokkcr-Planck model allow 
us to perform calculations in much wider ranges of parameter space, including varying the number of stars. 
This opens the possibility to compare Fokker-Planck calculations directly with direct iV-body calculations. 

With the improved Fokker-Planck model we have studied the evolution of globular clusters over a 
wide range of initial conditions. The surveyed initial parameters are identical to those adopted by CW. 
Our clusters live generally longer than the isotropic models of CW. The largest differences appear in cases 
where clusters dissolve promptly in the tidal field before core collapse occurs, except for cases of too prompt 
dissolution due to too severe mass loss via stellar evolution. For some initial conditions, our models reach 
core collapse, while the models of CW dissolve. The differences are rather small for long-lived clusters which 
start with steep mass functions and high concentrations and finally reach core collapse. These findings are 
consistent with the conclusions of Aarseth & Heggie (1998) who compare their iV-body models with the 
models of CW. 

About half of the clusters calculated in our survey are disrupted within 10 Gyr. Only four of our 36 
samples experience core collapse within 10 Gyr, and two of the four have lost more than 99% of their mass 
at the moment of core collapse. This suggests that the range of the initial conditions with which clusters 
reach core collapse within the age of the Universe is rather small in our survey. On the other hand, about 
20% of the known Galactic globular clusters are classified as post-collapse clusters (Djorgovski & King 1986; 
Chernoff & Djorgovski 1989; Trager et al. 1993). These post-collapse clusters might have been born with 
initial conditions that are preferable for core collapse at ~ 10 Gyr, e.g., a = 2.5 and 3 < W <^ 7. However 
it should be noted that not all of possible initial conditions of real globular clusters are covered by our (and 
CW's) survey. In particular the adopted initial mass functions are probably not the best choice. 

In the survey we find a few peculiar models which have very small masses (~ 1% of the initial mass) 
but are well concentrated at ~10Gyr. Such clusters happen to be classified as old open clusters, like 
Berkeley 17, rather than as globular clusters, depending on their luminosities, appearances, and positions 
relative to the Galactic disk. The present-day global mass functions of these model clusters are also peculiar 
in the sense that they decrease with decreasing stellar mass. A similar inverted mass function was observed 
in the globular cluster NGC 6712 by De Marchi et al. (1999). This indicates that NGC 6712 might have 
lost a very large fraction of its initial mass, if the observed mass function is really close to the global mass 
function and does not only reflect a peculiarity of the initial mass function of the cluster. 

Clusters which lose ^ 50% of their initial mass in a very short time scale may still survive for a 
long time. The principle that losing more than a half of the initial mass of a cluster inevitably results in 
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disruption is based on the virial theorem and on the assumption of impulsive mass loss (Hills 1980). For 
the initial conditions of W — 7, a — 1.5, and family 4, more than half the cluster mass is lost within a few 
half-mass crossing times (see Fig. ||e). But, in fact, still the cluster survives for more than ten billion years 
and finally experiences core collapse. 

Comparison between iV-body calculations with various implementations of the galactic tidal field and 
the cut-off boundary condition indicates that there may be a significant population of stars which are 
trapped outside the tidal radius of the cluster but are still bound to the stellar system. 
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Table 1. Stellar evolution models. 



TO ini /M 


log(fr 


Ms/yr) 


man 




cw 


SeBa 


CW 


SeBa 


0.40 


11.3 


11.3 


0.40 


0.40 


0.60 


10.7 


10.8 


0.49 


0.52 


0.80 


10.2 


10.4 


0.54 


0.56 


1.00 


9.89 


10.0 


0.58 


0.59 


2.00 


8.80 


9.02 


0.80 


0.69 


4.00 


7.95 


8.29 


1.24 


0.80 


8.00 


7.34 


7.64 


0.00 


1.33 


15.00 


6.93 


7.14 


1.40 


1.34 



Table 2. Families. 



Family 


-Few 


irix,cw/Gyr 


{{m 
Wo = 1 


>/M )(WGyr) 
Wo = 3 Wo = 


= 7 


1 


5.00 x 10 4 


128 


3.4 


2.5 


0.7 


2 


1.32 x 10 5 


339 


9.0 


6.5 


1.9 


3 


2.25 x 10 5 


577 


15.4 


11.0 


3.2 


4 


5.93 x 10 5 


1522 


40.6 


29.1 


8.3 



a The mean stellar mass (to) = 2.45M Q , 1.01M Q , and 0.66M Q , for 
a =1.5, 2.5, and 3.5, respectively. 
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Table 3. Survey results. 



VV o 










Family 














2 




3 










\-j VV 


A o 

a a 


o VV 


Aa 


CW 




o VV 


Aa 


1 


1.5 


D 


D 


D 


D 


D 


D 


D 


D 






9.2xl0~ 3 


0.017 


9.4xl0~ 3 


0.018 


9.3xl0~ 3 


0.018 


9.3xl0~ 3 


0.018 






0.78 


0.61 


0.78 


0.63 


0.78 


0.63 


0.79 


0.62 




2.5 


D 


D 


D 


D 


D 


D 


D 


D 






0.034 


0.31 


0.034 


0.32 


0.035 


0.32 


0.034 


0.33 






0.77 


0.55 


0.77 


0.56 


0.76 


0.56 


0.77 


0.57 




3.5 


D 


C 


D 


D 


D 


D 


D 


D 






2.5 


21.6 


2.9 


24.0 


3.1 


28.1 


3.2 


36.2 






0.64 


0.021 


0.73 


0.27 


0.74 


0.34 


0.76 


0.44 


3 


1.5 


D 


D 


D 


D 


D 


D 


D 


D 






0.014 


0.026 


0.014 


0.026 


0.014 


0.026 


0.014 


0.026 






0.53 


0.42 


0.55 


0.41 


0.55 


0.41 


0.54 


0.41 




2.5 


D 


D 


D 


D 


D 


D 


D 


D 






0.28 


2.2 


0.29 


2.7 


0.29 


2.9 


0.29 


3.0 






0.46 


0.27 


0.47 


0.32 


0.48 


0.34 


0.49 


0.37 




3.5 


C 


C 


C 


C 


D 


C 


D 


C 






21.5 


32.1 


44.4 


79.8 


42.3 


129 


43.5 


313 






0.078 


0.17 


0.035 


0.14 


0.085 


0.13 


0.28 


0.11 


7 


1.5 


D 


C 


D 


C 


D 


C 


D 


C 






1.0 


3.1 


3.0 


7.7 


4.2 


12 


5.9 


27 






0.022 


0.010 


3.3xl0~ 3 


0.012 


8.0xl0~ 3 


0.013 


0.023 


0.011 




2.5 


C 


c 


C 


C 


C 


C 


C 


C 






9.6 


10.1 


22.5 


23.8 


35.5 


37.7 


83.1 


88.7 






0.26 


0.32 


0.26 


0.32 


0.26 


0.32 


0.25 


0.31 




3.5 


C 


C 


C 


C 


C 


C 


C 


C 






10.5 


9.9 


31.1 


30.9 


51.3 


53.2 


131.3 


135 






0.57 


0.65 


0.51 


0.58 


0.48 


0.55 


0.49 


0.55 



a The results of Chernoff & Weinberg (1990, CW) are taken from their Tabic 5. Aa represents our 
resluts for the anisotropic Fokkcr-Planck models with the apocenter criterion. 

b The first entry describes the fate of the cluster at the end time of the simulation t cnd : C (core 
collapse) or D (disruption). The second entry is i cn a in units of 10 9 yr. The third entry is the cluster 
mass at £ en d, Af en d, in units of the initial mass. 
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Table 4. 


Compar: 


ison 


of CW, Ae, and Aa models for family 1. 


Wo a 


Model 


Fate i C nd 


M and 


t(M md ,cw) 


A/(*cnd.Cw) 








[Gyr] 


[Mo] 


[Gyr] 


[Mo] 


1 1.5 


cw 


D 


0.0092 


0.78 


0.0092 


0.78 




Ae 


D 


0.013 


0.69 


0.012 


0.98 




Aa 


D 


0.017 


0.61 


0.014 


0.98 


2.5 


CW 


D 


0.034 


0.77 


0.034 


0.77 




Ae 


D 


0.056 


0.61 


0.046 


0.85 




Aa 


D 


0.31 


0.55 


0.12 


0.89 


3.5 


CW 


D 


2.5 


0.64 


2.5 


0.64 




Ae 


D 


3.1 


0.45 


2.3 


0.60 




Aa 


C 


21.6 


0.021 


6.3 


0.82 


3 1.5 


CW 


D 


0.014 


0.53 


0.014 


0.53 




Ae 


D 


0.019 


0.49 


0.019 


0.76 




Aa 


D 


0.026 


0.42 


0.025 


0.82 


2.5 


CW 


D 


0.28 


0.46 


0.28 


0.46 




Ae 


D 


0.45 


0.38 


0.40 


0.58 




Aa 


D 


2.2 


0.27 


1.5 


0.72 


3.5 


CW 


D 


21.5 


0.078 


21.5 


0.078 




Ae 


C 


23.6 


0.097 




0.14 




Aa 


c 


32.1 


0.17 




0.37 


7 1.5 


CW 


D 


1.0 


0.022 


1.0 


0.022 




Ae 


D 


1.1 


0.030 


1.1 


0.046 




Aa 


C 


3.1 


0.010 


2.7 


0.12 


2.5 


CW 


c 


9.6 


0.26 


9.6 


0.26 




Ae 


c 


9.5 


0.27 








Aa 


c 


10.1 


0.32 




0.34 


3.5 


CW 


c 


10.5 


0.57 


10.5 


0.57 




Ae 


c 


9.7 


0.61 








Aa 


c 


9.9 


0.65 







a Most of the notation is defined in the notes to Tabic |3|. 

Ae represents anisotropic Fokkcr-Planck models with the energy criterion. 

c ^cnd,cw an d iend.cw denote the end mass and the end time, respectively, of each 
of the CW models. See text for more details. 
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Table 5. Comparison with TV-body models of FH and AH for family 1. 



Wo 


a 


Model 




tdis 




M dis 








[Gyr] 


[Gyr] 


[Mo] 


[Mo] 


1 


1.5 


FH3.7 




0.076 











Aa 




0.017 




0.61 




2.5 


FH4.0 




0.29 











Aa 




0.31 




0.55 




3.5 


FH4.1 




2.2 











Aa 


22 




0.021 




3 


1.5 


FH3.7 




0.11 











AH3.7 




0.075 











AH9.2 




0.12 











Aa 




0.026 




0.42 




2.5 


FH4.0 




0.91 











AH4.0 




3.0 











AH10.0 




3.4 











Aa 




2.2 




0.27 




3.5 


AH4.1 




>20 










AH10.4 




>20 










Aa 


32 




0.17 




7 


1.5 


AH3.7 




1.2 











AH9.2 




3.2 











Aa 


3.1 




0.010 






2.5 


AH4.0 


13 


>20 


0.21 








AH10.0 


11 


>20 


0.257 








Aa 


10 




0.32 






3.5 


AH4.1 


10.4 


>20 


0.606 








AH10.4 


10.5 


>20 


0.62 








Aa 


9.9 




0.65 





a Part of the notation is defined in the notes to Table 

b FH denotes iV-body models of Fukushigc & Hcggie 
(1995), and AH denotes TV-body models of Aarscth & Hcggie 
(1998) . A number associated with each of these model names 
indicates the galactoccntric distance R s (in kpc) at which 
that model is placed. 

c t cc and M cc denote the collapse time and mass, 
respectively, idis and Mdis denote the disruption time and 
mass, respectively. 
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Fig. 1. — Comparison of TV-body models (dashed lines) with Aa Fokker-Planck models with i/ csc —2, 2.5, 
and 3 (the lower, middle, and upper solid lines, respectively) for the evolution of the total mass. Larger f csc 
causes faster mass loss. Top panels: the initial conditions are (W , a, family) = (3, 2.5, 1). The average of 10 
runs is shown for the TV =1K TV-body model. Middle panels: (W , a, family) = (1,2.5,1). Bottom panels: 
(W , a, family) = (7,2.5,1). 



- 28 - 




Fig. 2. — Comparison between Aa Fokker-Planck models (solid lines) and TV-body models (dashed lines 
and dotted lines) for the evolution of the total mass. Thickness of the lines stand for largeness of 
N. The thickest solid line in each panel represents the Fokker-Planck model with the instantaneous 
escape condition, which corresponds to the the limit of TV —> oo. The other Fokker-Planck models are 
calculated with the crossing-time escape condition with ^ osc = 2.5 for finite TV. The dash-dotted lines 
represent the mass evolution expected when mass loss occurs only through stellar evolution (with no 
escapers). (a) The initial conditions are (W , a, family) = (3,2.5,1); TV =1K, 16K, and 32K, from right 
to left. Only for the TV =1K TV-body model, ±er/2 deviation from the mean of 10 TV-body runs is 
shown (two dotted-lines). (b) (W Q , a, family) = (1,2.5,1); TV=8K, 16K. (c) (W a , a, family) = (7,2.5,1); 
TV=4K, 8K. (d) (W , a, family) = (3,2.5,4); TV=8K, 16K. (e) (W , a, family) = (7,1.5,4); TV=8K, 16K. (f) 
(Wo, a, family) = (1,3.5,3); TV=8K, 16K. 
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Fig. 3. — Time evolution of the relative mass difference between Fokker-Planck and TV-body models, 
[Mp_p(t) — -^V-body W]/-^/V-body W ( son d lines), for the models shown in Fig. ||. The total mass of 
the TV-body models is also plotted (dashed lines). In each panel the initial conditions are denoted by three 
parameters (W Q , a, family). The thick lines represent the largest- N models in each panel. 
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Fig. 4. — Time evolution of the mass loss rates, £tot and £ csc (£tot > Cose), for the Fokker-Planck (solid lines) 
and TV-body (dashed lines) models shown in Fig. ^. In these plots time is expressed in units of the initial 
half-mass relaxation time, but for each panel the plotted period is the same as that in Fig. || (i.e., 3 Gyr for 
[b] and 10 Gyr for the others). Only the largest- iV models are shown for each set of (W , a, family), but in 
(a) N —IK models are also shown. 
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Fig. 5. — The global mass function for the model with the initial conditions of W = 7, a = 2.5, and family 1 
at t = 10 Gyr. The dashed line represents the initial mass function. The mass functions are normalized at 
each epoch. The Fokker-Planck results are given by symbols and the iV-body results are given by lines. The 
squares and the solid lines represent the mass functions for all the stars, and the circles and the dotted lines 
represent those for the main-sequence stars only. 
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Fig. 6. — Results of the Fokker-Planck survey simulations: Evolution of the total mass. Three panels from 
top to bottom show results for the initial mass functions of a=1.5, 2.5, and 3.5. Families 1 to 4 are plotted 
with solid, dotted, short dashed, and long dashed lines, respectively. The circles at the end points of some 
lines indicate that those simulations stop at the time of core collapse. 
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Fig. 7. — Results of the Fokker-Planck survey simulations: Evolution of the concentration c. The notation 
in this figure is the same as that in Fig. o. 
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a 

Fig. 8. — Graphical representation of the difference between the Aa models and the CW models in the 
survey results (see Table [|) for family 1 . The power-law index a of the initial mass function is given along 
the X-axis and the initial value for W is given along the Y-axis. The symbols C and D indicate the end 
states of the simulations, collapse and disruption, respectively. An X indicates that the result of CW is 
disruption where the Aa model reaches core collapse. The size of the letters is proportional to the lifetime 
at which we arrive as fraction of the lifetime of the cluster calculated by CW; a larger symbol indicates a 
larger discrepancy. 
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Fig. 9. — Comparison of the mass evolution between the Aa models of N — 3 x 10 5 clusters (dotted lines) 
and the Aa models of N — > oo clusters (solid lines, the models shown in Fig. |J). Only the models of family 
4 are shown. Three panels from top to bottom show results for the initial mass functions of a=1.5, 2.5, and 
3.5. 
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Fig. 10. — Mass as a function of time for various TV-body models with IK stars (for a — 2.5, W Q = 3, and 
family 1) for different implementations of the tidal field. Each line represents the mean of ten calculations. 
The solid line is calculated with a simple tidal cut-off as is used for all other TV-body calculations in this 
paper. The other lines give the results of calculations with a sef-consistent tidal field in which stars are 
removed beyond the radii r t (dotted), 2r t (short dashed), 10r t (long dashed), and 100r t (dash-dotted). The 
last two lines are almost indistinguishable. 
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Fig. 11. — The global mass function for the model with the initial conditions of W = 7, a = 1.5, and 
family 3, from the Fokker-Planck survey, at t = 12 Gyr. The open circles represent the mass function of the 
main-sequence stars, and the filled circles represent that of the white dwarfs. 



